######### This script will remove small patches from a species realized distribution ASCII
#########

#### Establish base directory

base.dir = '/home1/99/jc152199/MAXENT/'
setwd(base.dir)

#### Load library

library('SDMTools')

##### Need to loop through ASCIIs for multiple species
##### First identify ASCII directory

ascii.dir = '/home1/99/jc152199/MAXENT/output/'

##### Identify species

spp = dir(ascii.dir)

#### Remove directories which aren't skinks

spp = spp[c(1,13:15,17:19)]

###### Now loop through and identify individual ASCII for each species
##### Don't forget, this needs to be done for both the microCLIM, expCLIM, and BIOCLIM distributions

for (s in spp)

	{
	
	#### Species loop
	
	for (m in c('microCLIM','BIOCLIM'))
	
		{
		
		#### Model loop
		
		#### Read in an ASCII of realized distribution
		
		tasc = read.asc.gz(paste(ascii.dir,s,'/',m,'/output/',s,'_',m,'_Realized.asc.gz',sep=''))
		
		#### Keep an unchanged copy of tasc called base.asc
		
		base.asc =  tasc
		
		#### Convert to a binary matrix
		#### Don't need to apply MaxEnt thresh because that was done before producing the realized distribution
		
		tasc[which(tasc>0)]=1
		
		#### Run class-stat on tasc to determine total area for the species distribution
		
		cstat = ClassStat(tasc,cellsize=250,bkgd=NA,latlon=TRUE)[2,]
		
		### Extract the total area of the species realized distribution
		
		s.area = cstat[1,3]
		
		#### Run ConnCompLabel to identify individual patches
		
		pasc = ConnCompLabel(tasc)
		
		#### Run Patch Stat on pasc
		
		pstat = PatchStat(pasc, cellsize=250, latlon=TRUE)
		
		### Calculate the area of each patch as a proportion of the total possible area
		
		pstat$propreal = (pstat$area/s.area)*100
		
		#### Identify patches less than 5% of total area
		
		pstat = pstat[which(pstat$propreal>=.05),]
		
		#### Now use a loop to find all positions in pasc that match a patchID in pstat
		
		pos2write = NULL
		
		for (p in pstat$patchID[which(pstat$patchID>0)])
		
			{

			tpos = as.data.frame(which(pasc==p, arr.ind = T))
			
			pos2write = rbind(tpos,pos2write)
			
			}
		
		### Extract data from positions in pos2write, first get lat/longs to extract from
		
		pos2write$lat = getXYcoords(tasc)$x[pos2write$row]
		pos2write$long = getXYcoords(tasc)$y[pos2write$col]
		
		#### Extract data from the base ASCII
		
		pos2write$data2write = extract.data(cbind(pos2write$lat,pos2write$long),base.asc)
		
		#### Convert finite values in base.asc to zero
		
		base.asc[which(is.na(base.asc)==F)]=0
		
		#### Write data from pos2write back onto base.asc
		
		base.asc[cbind(pos2write$row,pos2write$col)]=pos2write$data2write
		
		#### Write out base.asc
		
		write.asc.gz(base.asc,file=paste(ascii.dir,s,'/',m,'/output/',s,'_',m,'_Realized_No_Small_Patches',sep=''))
		
		}
		
	cat('\n',s,' - Completed\n',sep='')	
		
	}
	
#### Done

#### Checking

realized = read.asc.gz(paste(ascii.dir,s,'/',m,'/output/',s,'_',m,'_Realized.asc.gz',sep=''))
realizednopatches = read.asc.gz(paste(ascii.dir,s,'/',m,'/output/',s,'_',m,'_Realized_No_Small_Patches.asc.gz',sep=''))

realized[which(realized>0)]=1
realizednopatches[which(realizednopatches>0)]=1

rccl = ConnCompLabel(realized)
rnpccl = ConnCompLabel(realizednopatches)
		
rstat = PatchStat(rccl, cellsize=250, latlon=TRUE)
rnpstat = PatchStat(rnpccl, cellsize=250, latlon=TRUE)

### Looks like it works
### realized has 457 patches
### realizednopatches has only 25
		
		

		



